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Abstract 

We investigate the possibility that freely rotating cylinders with an aspect 
ratio L/D = 0.9 exhibit a cubatic phase similar to the one found for a system 
of cut-spheres. We present theoretical arguments why a cubatic phase might 
occur in this particular system. Monte Carlo simulations do not confirm 
the existence of a cubatic phase for cylinders. However, they do reveal an 
unexpected phase behavior between the isotropic and crystalline phase. 

I. INTRODUCTION 

The phase behavior of apparently simple systems can be surprisingly rich. In 1957 Alder 
and Wainwright reported the first simulations of a system of hard spheresB. They showed 
that, upon increasing density, this system exhibits a phase transition from an homogeneous 
liquid to a crystalline FCC phase, even though there is no attractive interaction. 

An almost equally simple system that exhibits an even richerphase behavior, is the hard 
spherocylinder model, that was introduced by Onsager in 19491. Onsager used this model 
to show that hard, elongated, rodlike particles must undergo a transition from the isotropic 
to nematic phase. Subsequently, computer simulations, revealed that a smectic A and two 
different crystalline phases can existFi, depending on the aspect ratio. Recent extensions 
of the Onsager theory to higher densities and smaller aspect ratios can account for these 
novel phasesOia. 

In order to study the behavior of disklike objects Eppenga and Frenkel performed simu- 
lations of a system of infinitely thin disks0. By construction, the packing fraction of such a 
system is zero. As a consequence, after aligning the particles, the system can easily be com- 
pressed, nevertheless, the nematic to columnar transition for this system has been studiecfi 
A different model for disklike object was proposed by Veerman and Frenkefi They used 
cut-spheres, which are obtained by symmetrical slicing off the polar caps of a sphere. As a 
function of the length over diameter ratio of these particles, they connect the infinitely thin 
discs with the spheres. Surprisingly, cut-spheres with an aspect ratio of about 0.2, appear to 
exhibit an orient at ionally ordered phase with a cubic symmetry, which is called the cubatic 
phase, even though these particles themselves have a cylindrical symmetry. 

In this article we investigate a possible explanation of the existence of the cubatic phase 
for cut-spheres. To this end we study, by computer simulation, a system of freely rotating 
cylinders. In Sec. |TJ we present a theoretical basis for looking at cylinders. In Sec. |U 
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we elucidate some of the techniques used in our Monte Carlo simulation. This involves 
the overlap criterion of these particles and an additional type of move that we used in 
our simulations. In Sec. [IV] we show the main results from the computer simulations and 
discuss the unexpected behavior of cylinders between the isotropic and crystalline phase. 
The symmetry of this new phase is analyzed in Sec. [VJ In Sec. |VT] we determine the 
coexistence between the isotropic liquid and the crystalline phase by means of a free energy 
calculation. In Sec. |V11| we discuss the main results and give some suggestions for future 
research. 



II. CUBATIC PHASE 

The cubatic phase is a long-range orient at ionally ordered phase without any positional 
order of the particles. In a uniaxial nematic phase the particles, e.g. rods, tend to align along 
a single preferred direction. In a biaxial phase for single component systems, e.g. biaxial 
ellipsoidstj, the different molecular axes of the particles align in three different orthogonal 
directions. In the special case that these three directions are also equivalent, in the sense 
that each molecular axis is with equal probability in any of the three directions, we obtain 
an orientationally ordered phase with cubic symmetry. In the absence of translational order 
this is the elusive cubatic phase. 

Evidence for this cubatic phase was first found in computer simulations of cut-spheres 
with an aspect ratio L/D = 0.20. However, subsequent theoretical work indicates that 
particles, consisting of three perpendicular, elongated rods of approximately the same length, 
at least in the Onsager limit of infinite aspect ratios, should also form this cubatic phase@'0. 
Experimentally, however, this cubatic phase has never been observed. 

Figure [1| shows a snapshot of the cubatic phase of cut-spheres. The figure suggests that 
the disklike objects tend to form stacks of several particles in an approximately cylindrical 
shape. These cylindrical clusters are then arranged in such a way that there is an overall 
cubatic order. As a first attempt to explain the cubatic phase for cut-spheres we will attempt 
to describe the system on the level of these aggregates. This is a simplification, because the 
stacks in Fig. [I] are obviously not perfect cylinders but polydisperse in both length and 
shape. Yet, in what follows, we model this system as a collection of monodisperse, hard 
cylinders. 

In order to select the optimal aspect ratio that might give rise to a cubatic phase, we 
first consider the excluded volume of two cylinders. Onsager already gave an expression for 
the excluded volume S of two arbitrary cylinders with lengths Lj and diameters D$ 

£(7) = \D 1 D 2 {D 1 + D 2 ) sin 7 + L X L 2 (D X + D 2 ) sin 7 
+ L x (f D\ + D 1 D 2 E{sm 1 ) + f D\\ cos 7 |) 

+ L 2 (I D\ + D 1 D i E(sm 7 ) + \ D\\ cos 7 |) (1) 

where 7 is the angle between the two main axes of the cylinders and E(x) is the complete 
elliptical integral of the second kind. 

In Fig. we have plotted the excluded volume of two identical cylinders as function 
of the angle 7, for several aspect ratios. To facilitate comparison we normalized the ex- 
cluded volume to unity for perpendicular orientations. Note that for the aspect ratios, zero 
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and infinity, £ (7) reduces to the same monotonically increasing function. This reflects the 
tendency of strongly anisotropic particles to align. For an aspect ratio of the order unity, 
however, we see that the maximum excluded volume is found between the parallel and per- 
pendicular orientation. In fact, for the latter two orientations the exclude volume exhibits 
a minimum, although the deepest minimum always corresponds to a parallel orientation. 

In a cubatic phase both parallel and perpendicular orientations should be present. How- 
ever, it is logical to look for that aspect ratio that minimizes the relative difference between 
both minima. The fact that parallel orientations lead to a smaller excluded volume, and 
are therefore favored, could be compensated by the fact that forming a cubatic phase would 
lead to a higher orientational entropy. This is because particles in a cubatic phase can be 
aligned in three different directions instead of a single one as in the case of only parallel 
orientations. In Fig. |3] we have plotted the ratio of the perpendicular and parallel orien- 
tation as a function of the aspect ratio. The minimum of the curve occurs for the aspect 
ratio L/D = ~ 0.886, where the excluded volume of perpendicular orientations is 

only 1.133 times larger than that of parallel orientations. 

III. SIMULATION 

In our search for a possible cubatic phase for hard cylinders, we performed Monte Carlo 
simulations of a system of freely rotating, hard cylinders with an aspect ratio L/D = 0.9. 

A. Overlap criterion 

an overview is given of the most widely used overlap criteria for convex hard particles. 
The one for cylinders is only outlined. It consists out of three steps. 

1. spherocylinder overlap: For the first step we replace the cylinders by spherocylinders 
by adding hemispheres on both sides of the cylinders. We then compute the closest 
distance between the spherocylinders. If the spherocylinders do not overlap the cylin- 
ders will not either. If the spherocylinders do overlap we can only conclude that the 
cylinders overlap if the vector of shortest distance intersects both cylindrical parts. 

2. disk-disk overlap: If the last condition is not satisfied we need to continue with the 
second step. In this step we check whether there is an overlap of the flat faces of the 
cylinders. 

3. disk-cylinder overlap: If no overlap is found in the second step, and the spherocylinders 
in the first step did overlap we need to proceed with the third step, which is the most 
time-consuming one. We have to check whether one of the flat faces of a cylinder 
overlaps with the cylindrical hull of the other cylinder. 

According toEil this last problem reduces to calculating overlaps between two planar 
ellipses, and hence is a special case of the overlap between two ellipsoids. It is this last step 
which is the most difficult. 

For ellipsoids there exist two different criteria to determine whether there is an overlap 
or not. The first is due to Vieillard-Baron0 and the second due to Perram and WertheimB. 
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Both criteria determine whether there exists an overlap without calculating any points which 
both particles have in common. But although the intersection of the plane in which the flat 
face lies with the cylindrical hull of the other particle is an ellipse, the ellipse is not necessary 
complete. This happens if the plane also intersects with one of the flat faces of the second 
cylinder. In that case we need to check whether there is an overlap with the incomplete 
ellipse. 

The alternative route is to look at the problem as the overlap between a flat face of one 
cylinder with a sequence of circles of the second particle. This results in a fourth order 
polynomial equation and can in principle be solved analytically. The main problem with 
this approach is that during simulations it will result in severe problems, and not always 
find the correct outcome. The reason for this is that the coefficients of the polynomial 
equation are inversely proportional to the sine of the angle between the particles. Therefore 
the coefficients can differ by orders of magnitudes, which makes it hard to solve and not 
always accurate. 

The method we applied for the last step in the overlap criterion therefore is a direct 
minimization of the distance between a point on the edge of the first cylinder and a point 
on the axis of the second cylinder. Although it is slower, it has proven to be very reliable. 

B. Flip-move 

In Monte Carlo simulations we wish to sample the complete configuration-space in an 
efficient way. In our Monte Carlo simulations we keep the number of particles N, the pressure 
P and the temperature T are fixed. We performed several distinct type of trial moves. We 
try to translate and rotate particles by a small amount. In addition we perform volume 
changing moves in order to equilibrate the density under the applied pressure. For details 
about these standard techniques, we refer to0. 

At high densities, we can only perform relatively small trial moves, because otherwise 
the acceptance of the moves becomes negligibly small. However, this may create a kinetic 
barrier to equilibrate a cubatic phase, in which particles are parallel or perpendicular to each 
other, while intermediate orientations are less likely due to steric hindering. It is difficult to 
achieve a rotation of a cylinder over 90° by a succession of small rotational moves. 

In order to sample the configuration space of this system more efficiently, we introduce an 
extra trial move: the flip move. Rather than waiting for the rare reorientation of a cylinder 
over 90° through a succession of small angular jumps, we attempt to rotate it directly to 
a perpendicular orientation. The flip-move is a proper Monte Carlo move that satisfies 
detailed balance. By construction it is such, that it does not break the symmetry of an 
isotropic phase. In the case of short cylinders this trial move is surprisingly effective. 

To illustrate this we have shown in Fig. |] three steps of the continuous rotation of a 
short cylinder in a dense ordered phase. The first picture shows that for parallel particles 
there is still some space to move. However, if the cylinder in the middle is rotated by 45°, 
its dimension in the plane changes causing it to overlap with neighboring particles. As the 
rotation continues the excluded volume of the cylinder decreases again. By including the 
flip-move we bypass the unfavored intermediate state. 

During our simulations on average 1 out of every 10 rotations will be an attempt to flip 
a particle. To this end, we select at random an orientation in the plane perpendicular to 
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the direction of the particle and accept it if there is no overlap in the new situation. The 
acceptance of this move depends of course strongly on the aspect ratio and density, but can 
be as large as 5%, even in a crystalline phase. 



C. Nature of crystalline phase 

There are several relevant crystalline structures to study for a system of cylinders with an 
aspect ratio of order unity. The highest possible density in a system of cylinders is reached 
when all cylinders are perfectly aligned in a close-packed structure with the particles ordered 
in hexagonal layers. This leads to a maximum packing fraction <fi = n/\/l2. 

In the case of spheres there are many different possible ways to stack the layers to form 
a regular crystal. The simplest stackings are the face centered cubic (FCC) and hexagonal 
close-packed (HCP) structure. In both structures, the close-packed layers are shifted with 
respect of each other, such that the centers of mass in one layer are above the holes in the 
layer below. This leads to three different layer positions labeled by A, B and C. The FCC 
crystal corresponds to ABC stacking, while ABAB stacking is characteristic for the HCP 
crystal^. 

In the case of a densely packed crystal of cylinders, there is an infinity of possible stackings 
positions, provided that the close-packed planes are flat. However away from close packing, 
the system will prefer a situation where the layers are stacked in an AAA fashion. Other 
stackings are not stable and will relax to this structure. 

If the cylinders are not confined to close-packed planes, the AAA crystal phase could 
deform in a columnar phase. In this phase the cylinders are organized in a hexagonal array of 
columns. However, there is no long-range correlation between the longitudinal displacement 
of different columns. 

Another crystalline phase that is possibly relevant is the simple cubic crystal, in which the 
orientation distribution of the cylinders has the same symmetry as the lattice. Although this 
phase cannot be close packed, it allows us prepare a structure that has the same orientational 
symmetry as a cubatic phase. This is achieved by aligning the cylinders at random with any 
of the three axes of the crystal. 



IV. SIMULATION RESULTS 

Constant NPT-MC simulations were performed on a system of 720 particles. Long 
runs were needed to collect good statistics (typically 10 5 trial moves per particle excluding 
equilibration). On average a volume change was attempted N times less frequently than the 
simple particle moves. One in every ten trial rotations was a flip-move. 

In our NPT simulations of the solid phase, we allowed the shape of the simulation box 
to change to an arbitrary parallelepiped rather than keeping it rectangular. Such a move 
was first used by Parrinello and Rahmanihi in MD simulations and later by Najafabadi 
and Yip in MC simulations". Allowing for shape changes is only useful in crystalline 
structures, because then the presence of the crystal structure acts as a restoring force and 
will ensure that the box shape cannot become extremely anisotropic. In a liquid phase, 
extreme deformations could occur, because a natural restoring force does not exists. 
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The resulting equation of state is shown in Fig. ||, where we plotted the reduced pressure 
(5Pv versus the packing fraction 0, where v is the volume of a single cylinder and j3 = 
1/(/jbT) the inverse temperature. The two branches are obtained by compressing a low- 
density liquid phase and expanding a high-density AAA-crystal. The equation of state 
provides no indication that a cubatic phase, or any other intermediate phase, exists. 

If we start with a crystalline phase in which the hexagonal layers are stacked in a different 
fashion, for instance ABC, the layers will swiftly slide with respect to each other to form an 
AAA-stacking. Occasionally, we find that the columns shift with respect to each other, but 
this is probably a finite size effect, rather than an indication of a columnar phase. 

We can also start from a simple cubic crystal in which the particles are either regularly or 
randomly aligned with respect to the three axes of the crystal. This phase is also not stable 
but tends to transform into an the AAA-crystal. It is observed that this sometimes leads to 
two AAA-crystals with different orientations. However, it is likely that, given enough time, 
the grain boundary that has formed will anneal out. 

The flip- move is extremely effective. Even for a packing fraction = 0.7 one out of 
every 1000 attempted flips is accepted. For = 0.65 this is already one out of every 
100 and increases to about 7% near the transition to the isotropic phase. A simulation, 
without using this flip-move, showed that at = 0.65 particles can achieve a flip via a 
continuous rotation, but this occurs with a very low probability, only once in a simulation 
of 10 5 sweeps. Therefore the introduction of the flip-move speeds up this slow 'dynamical' 
process enormously, leading to a faster convergence to equilibrium. 

In order to get an impression what these systems look like, three snapshots are shown 
in Fig. ^. The first snapshot shows a high-density isotropic phase. The second represents 
a high-density, almost perfect AAA-crystal at a packing fraction = 0.684, with one orien- 
tational defect. The third and most interesting snapshot is obtained at the low-density side 
of the crystalline branch. This is the place where, the cubatic phase would be observed if it 
existed. The phase that we observe appears to be neither isotropic nor perfectly crystalline. 
A stunning feature of this phase is the fact that the shape of the simulation box changes 
from a rectangular shape at high densities to a non-orthogonal shape in the intermediate 
region. As the symmetry of the structure differs from both the high density solid and the 
low density liquid, we must classify it as a separate phase. Since the equation of state shows 
no hint of a first-order phase transition, the novel phase is probably joined to the crystalline 
phase by either a continuous or weakly first-order transition and is found at packing fractions 
between 0.55 and ~ 0.65. 



In order to analyze the nature of the intermediate phase we discuss the translational 
and orientational order separately. To characterize the positional order we use the bond 
orientational order parameters Qi used by Steinhardt et aZ.il, 



V. SYMMETRY OF THE INTERMEDIATE PHASE 
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where we determine the polar angles of the vector r connecting neighboring particles. Only 
particles within a distance of twice the length of the particles are used. 

The bond orientational order parameters for / = 4, 6 and 8 are shown in Fig. [F[ As 
can be seen from the figure, the behavior of these order parameters does not indicate any 
sudden structural phase transition. Figure [8] shows the projections of the centers of mass 
of the particles close to melting (<fi = 0.567). The top view shows clearly the hexagonal 
structure, indicating the presence of columns. The side view shows that columns lie within 
planes. 

We use the system at packing fraction <fr = 0.567 to illustrate the orientational order 
in the phase near the transition. In Fig. |9] we plot all orientations of the particles on a 
unit sphere. It is clear that particles are oriented either along a main axis or in a plane 
perpendicular to that direction. 

In order to detect the orientational order we calculate the following matrix 

1 N 

Q = a7 E - I 1 ) > ( 3 ) 

i=l 

where Ui is the director of a particle and N is the total number of particles. The nematic 
order parameter P2 is defined as the average value of the maximum eigenvalue of the matrix 
Q and the nematic direction as the corresponding eigenvector. Furthermore we evaluate the 
order parameters //, defined by 




^a, m qj ). (4) 



These order parameters are defined in terms of the modified spherical harmonics Cj, m , and 
are frame independent. I2 is related to the nematic order parameter. But whereas the 
nematic order parameter P2 measures uniaxial order only, I2 measures the combination of 
uniaxial and biaxial order. The order parameter 7 4 is also able to measure uniaxial order, 
but now via P 4 , and biaxial order. In addition, however, it is also able to measure cubatic 
order. It is therefore a non-zero I4 combined with the absence of nematic order, and thus a 
zero valued I2, which would indicate the presence of cubatic order in the system. 



Both order parameters, P2 and I4, are shown in Fig. [10| The non-zero value of P2 on the 
complete crystalline branch, indicates that there is a single preferred or average direction 
in the system. It turns out that this direction is along the direction of the columns of the 
crystal. For low density crystals, however, the value of P2 is small and comparable to the 
nematic order close to the isotropic-nematic phase transition. This is explained by the fact 
that particles are either aligned along the columns or perpendicular to them. For lower 
densities more particles will have there direction perpendicular to the columns. It is not 
favorable to have a direction in between the two extreme orientations. 

In order to detect the orientational order in the plane perpendicular to the columns, we 
finally consider the average of the functions sin Z 93, where the angles ip of the directors of the 
particles are measured in the plane perpendicular to the average direction obtained from (||) 

Ox = (|C M |>- (5) 
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These order parameters are also shown in Fig. |TTj, and indicate that there is a small, 4-fold 
symmetry present. This is caused by the preference of particles in a layer to be either parallel 
or perpendicular. The 4-fold symmetry is, however, in conflict with the 6-fold symmetry 
of the crystal. These conflicting symmetries are probably the reason of the non-orthogonal 
unit cell. 



In Fig. |ll| we show the radial distribution function g(r) and the orientational correlation 



functions gi(r) for I = 2 and 4, that are defined by 

9l (r) = (P,(u(0).u(r))), (6) 

where Pi is the Ith Legendre polynomial and u(r) is a unit vector along the axis of a 
particle at distance r from the reference particle. These correlation functions are a measure 
of the range of the orientational order. As can be seen from Fig |ll|, g 2 (r) and g±{r) are 
finite throughout the complete simulation box. Although the correlation function g±{r) 
is measuring a combination of nematic-like and cubatic-like order, the magnitude of g±{r) 
compared to that of #2(7) indicates that the cubatic-like order is dominant. This is different 
from a normal nematic phase where g^ (r) is larger that g±{r). 

Finally we show in Fig. [T2| the density correlation functions h a (r) defined by 

hair) = (p(0)p(r))/p 2 , (7) 

where the distance r is measured along the a direction and in units of the simulation box 
length. The functions h x and h y show distinct features, as could be expected from the 
projections shown in Fig. |8|, and the number of peaks correspond to the number of layers 
present in that direction. In the z direction a weak modulation is observed, although the 
peaks of the function h z are less pronounced. This can partly be understood if one realizes 
that the fluctuations in the positions along the columns can be of the order of 20% of the 
effective size of the particles in that direction. These fluctuations are even enhanced by the 
fact that particles are oriented along and perpendicular to the columns, but the effective 
size in both cases is different. For the system sizes that we studied, the density modulation 
in the z direction persists over the entire simulation box. However, we cannot rule out the 
possibility that the modulation will decay in even larger systems. If this were the case, the 
system would, in fact, be in a columnar phase. 

In summary, in our simulations we find evidence for an intermediate phase with 3D- 
crystalline structure in which hexagonal layers are stacked in an AAA fashion. The orien- 
tations of the particles are coupled to the crystal axes and are either along the columns or 
within the layers, in which case they have a slight 4-fold symmetry. Neither the equation of 
state, nor the order parameters show any indication that there is a first order phase tran- 
sition along the crystal branch to the perfect, aligned crystal. The novel phase is found at 
packing fractions between m 0.55 and ~ 0.65, although this upper boundary is some- 
what arbitrary as we found no clear phase transition. Above this boundary, however, less 
than 1% of the particles is misaligned with the crystal. 



VI. FREE ENERGY CALCULATION 

In order to determine the coexistence between the isotropic liquid and the ordered phase 
we have to find points on the equation of state with equal pressure and chemical potential. 
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To this end we need to determine the free energy. By performing a thermodynamic inte- 
gration along the equation of state we can evaluate the free energies up to a constant. This 
integration is only allowed if there is no first order transition along the integration path, 
which in our case is true because all evidence indicates that the transition on the crystal 
branch is continuous. 

In order to proceed we need to determine one reference point on either branch of the 
equation of state. In the isotropic phase this is easy because we can use the ideal gas as a 
reference state 

F{p) = FM + f dp' P{p,) ~ P ' . (8) 
Jo P 

For the crystal phase we use as a reference system an Einstein crystal with the same 
structure^, which in our case is a perfectly aligned AAA-crystal. We connect the two 
systems by a one parameter Hamiltonian which consists of two parts, the coupling of the 
particles to their equilibrium lattice positions and the alignment for the orientations 

Hx = A J> - r?) 2 + A si* 2 ^), (9) 

i i 

where A is the coupling parameter, ff are the lattice positions, ri the positions of the particles 
and 9i are respectively the orientation of the particles with respect to the preferred direction. 
By slowly increasing the value of A the system will order according to the imposed field. We 
used here the same coupling parameter for both fields, but one can also use different values 
to do the orientational and positional ordering separately. 

The free system with A = can be related to the Einstein crystal for which A ^> 1 by 

(3F(p) (3F em {\ max ) f Xmax r X ^ X 2^ f Xmax ^^-2n^ ^gV 

' dX < or >\ — d\ < sin > A — — , (10) 



N N J J N 

where < 5r 2 >\ is the mean square displacement and < sin 2 6 >\ the average sine squared 
of the angle between the directors of the particles and the direction of the alignment field. 
The last term corrects for the fact that the center of mass during this simulation needs to be 
fixed. The value of X max is chosen such that a system using this Hamiltonian only will not 
cause any overlaps in the system. For that limit we can derive the value of the free energy 
of the Einstein crystal@ 

(3F em {X) = -\ logiV - \{N - 1) log(— ) - N\og(—). (11) 

In the case that overlaps do occur, it is possible to correct for thisS. By simulating the 
system for different values of A the integrals in (|10"D can be evaluated numerically. In order 
to minimize the error we used the Gauss- Legendre quadrature. 

The coexistence is obtained is at packing fraction = 0.515 ± 0.003 for the isotropic and 
= 0.552 ± 0.006 for the crystalline structure, at a pressure {3Pv = 8.76 ± 0.08 (Fig. ||). 
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VII. DISCUSSION 



Our study of the hard cylinder system was inspired by the possibility that it might exhibit 
a cubatic phase. We have presented Monte Carlo simulation for cylinders with an aspect 
ratio L/D = 0.9. In the standard NPT-simulations we introduced the flip-move, which 
rotates a particle directly to a perpendicular orientation. Furthermore we allowed the box 
shape to change to an arbitrary parallelepiped to facilitate equilibration of the crystalline 
phase. 

At high pressures the system is found in an AAA-crystal with particles aligned along the 
direction perpendicular to the layers. By lowering the pressure, particles reorient, deform 
the perfect crystal and form an intermediate phase in which the 6-fold symmetry of the 
crystal is slightly distorted and only very weak layering can be observed. The orientations 
are either along the columns or perpendicular to them. Towards the phase transition to the 
isotropic liquid, a weak, but significant, 4- fold orientational order in plane develops. There 
is no indication that this structural phase transition is first order. The phase coexistence 
between the isotropic and crystal phase was determined at packing fractions <fi = 0.515 in 
the isotropic phase and <fi = 0.552 for the crystal phase. 

Our simulations did not reveal a true cubatic phase. In fact, a theoretical analysis of this 
systemEi suggests that the isotropic to cubatic phase transition occurs at a density beyond 
that of the isotropic to nematic phase transition. This theory predicts that the isotropic 
to cubatic transition would take place at a packing fraction rs 0.80. However, at this 
density positional order cannot be neglected and should be incorporated in the theoretical 
description. 

The present results suggest that it is an oversimplification to represent stacks of cut- 
spheres by monodisperse hard cylinders. As can be seen in Fig. [I], the stacks are polydisperse 
in length and shape. It might be possible that a cubatic phase is stabilized by introducing 
polydispersity to the system. However, such a study is outside the scope of the present 
paper. 

ACKNOWLEDGMENTS 

We thank Martin Bates for a critical reading of the manuscript. The work of the FOM 
Institute is part of the research program of FOM and is made possible by financial support 
from the Netherlands Organization for Scientific Research (NWO). 



10 



REFERENCES 



1 B. J. Alder and T. E. Wainwright, J. Chem. Phys. 27, 1208 (1957). 
2 L. Onsager, Ann. (N.Y.) Acad. Sci. 51, 627 (1949). 

3 A. Stroobants, H. N. W. Lekkerkerker, and D. Frenkel, Phys. Rev. Lett. 57, 1452 (1986). 

4 P. Bolhuis and D. Frenkel, J. Chem. Phys. 106, 666 (1997). 

5 D. Frenkel, H. N. W. Lekkerkerker, and A. Stroobants, Nature 332, 822 (1988). 

6 A. M. Somoza and P. Tarazona, Phys. Rev. A 41, 965 (1990). 

7 A. Poniewierski and R. Holyst, Phys. Rev. A 41, 6871 (1990). 

8 H. Graf, H. Ldwen, and M. Schmidt, Prog. Colloid Polym. Sci. 104, 177 (1997). 
9 H. Graf and H. L6wen, Phys. Rev. E 59, 1932 (1999). 
10 R. Eppenga and D. Frenkel, Mol. Phys. 52, 1303 (1984). 

11 M. A. Bates and D. Frenkel, Phys. Rev. E 57, 4824 (1998). 

12 J. A. C. Veerman and D. Frenkel, Phys. Rev. A 45, 5632 (1992). 

13 P. J. Camp and M. P. Allen, J. Chem. Phys. 106, 6681 (1997). 

14 D. Frenkel, in Liquids, Freezing and Glass Transition, edited by J. P. Hansen, D. Levesque, 
and J. Zinn- Justin (North-Holland, Amsterdam, 1991), pp. 689-762. 

15 R. Blaak and B. M. Mulder, Phys. Rev. E 58, 5873 (1998). 

16 M. P. Allen, G. T. Evans, D. Frenkel, and B. M. Mulder, Adv. Chem. Phys. 86, 1 (1993). 

17 J. Vieillard-Baron, J. Chem. Phys. 56, 4729 (1972). 

18 J. W. Perram and M. S. Wertheim, J. Comp. Phys. 58, 409 (1985). 

19 D. Frenkel and B. Smit, Understanding Molecular Simulation. From Algorithms to Appli- 
cations (Academic Press, Boston, 1996). 

20 C. Kittel, Introduction to Solid State Physics, 6th ed. (John Wiley & Sons, New York, 
1986). 

21 M. Parrinello and A. Rahman, Phys. Rev. Lett. 45, 1196 (1980). 

22 M. Parrinello and A. Rahman, J. Appl. Phys. 52, 7182 (1981). 

23 M. Parrinello and A. Rahman, J. Chem. Phys. 76, 2662 (1982). 
24 R. Najafabadi and S. Yip, Scripta Metall. 17, 1199 (1983). 

25 P. J. Steinhardt, D. R. Nelson, and M. Ronchetti, Phys. Rev. B 28, 784 (1983). 

26 D. Frenkel and A. J. C. Ladd, J. Chem. Phys. 81, 3188 (1984). 

27 J. M. Poison, E. Trizac, and D. Frenkel, (Unpublished). 

28 R. Blaak, Ph.D. thesis, University of Utrecht, 1997. 



11 



FIGURE CAPTIONS 



1. Snapshot of the cubatic phase for cut-spheres. 

2. The excluded volume of two identical cylinders as function of the mutual angle 7 
for different aspect ratios L/D = and 00 (solid), 1 (dashed), 2 (long dashed) and 
10 (dotted-dashed). The excluded volume is normalized to unity for perpendicular 
orientations. 

3. The ratio of excluded volume for perpendicular and parallel orientations as a function 
of the aspect ratio of the cylinders. The minimum is found for L/D = yfjr/2. 

4. Top view of a hexagonal structure for cylinders with aspect ratio of order unity. The 
left figure shows the aligned case and the middle figure the case when the particle is 
perpendicular to a layer. The right figure shows an intermediate orientation causing 
overlap. 

5. The equation of state of cylinders with L/D = 0.9 for the isotropic (o) and crystal 
(o) phase. The reduced pressure (3PVq, where vo is the volume of a single particle, is 
plotted versus the packing fraction 0. The point of coexistence between the isotropic 
and crystal branch are connected by the solid line. 

6. Three snap shots from simulations, corresponding to isotropic (top left), high den- 
sity crystal (top right) and ordered phase near the transition (bottom). The packing 
fractions for these phases are respectively 0.478, 0.684 and 0.567. 

7. The bond orientational order parameters Q^(o), Qq(o) and Qs(a) as function of the 
packing fraction in the ordered phase. 

8. The projections of a configuration at packing fraction = 0.567. The top view (top) 
shows the hexagonal pattern of the columns which can be seen from the side (under 
left). In the front view (right under) planes are more difficult to observe. 

9. The orientational distribution for a configuration at = 0.567. 

10. The orientational order parameters as function of the packing fraction in the ordered 
phase. The nematic (o), cubatic (o) and C 4 (□) order parameters are non-zero, C 2 (a) 
and C 6 (v) are negligible small. 

11. The radial distribution function g(r) (o) and the orientational correlation functions 
92(f) (□) and g±(r) (o) as a function of the distance r measured in units of the diameter 
of the cylinders at = 0.567. 

12. The correlation functions h a (r) as function of the distance along the axis and measured 
in units of the box length at = 0.567. h x (solid) and h y (dotted-dashed) clearly 
indicate the presence of layers in the x respectively y direction. Although h z (dashed) 
has less pronounced peaks, also here layering is visible. 
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FIGURES 




FIG. 1. R. Blaak, Journal of Chemical Physics 
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FIG. 2. R. Blaak, Journal of Chemical Physics 



14 




FIG. 3. R. Blaak, Journal of Chemical Physics 
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FIG. 4. R. Blaak, Journal of Chemical Physics 
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FIG. 5. R. Blaak, Journal of Chemical Physics 
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FIG. 6. R. Blaak, Journal of Chemical Physics 
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FIG. 7. R. Blaak, Journal of Chemical Physics 
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FIG. 8. R. Blaak, Journal of Chemical Physics 
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FIG. 9. R. Blaak, Journal of Chemical Physics 
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FIG. 10. R. Blaak, Journal of Chemical Physics 
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FIG. 12. R. Blaak, Journal of Chemical Physics 
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